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Abstract 

Based on a non-equilibrium mechanism for spatial pattern formation we study 
how position information can be controlled by locally coupled discrete dynamical 
networks, similar to gene regulation networks of cells in a developing multicel- 
lular organism. As an example we study the developmental problems of domain 
formation and proportion regulation in the presence of noise, as well as in the 
presence of cell flow. We find that networks that solve this task exhibit a hi- 
erarchical structure of information processing and are of similar complexity as 
developmental circuits of living cells. Proportion regulation is scalable with sys- 
tem size and leads to sharp, precisely localized boundaries of gene expression 
domains, even for large numbers of cells. A detailed analysis of noise-induced 
dynamics, using a mean-field approximation, shows that noise in gene expression 
states stabilizes (rather than disrupts) the spatial pattern in the presence of cell 
movements, both for stationary as well as growing systems. Finally, we discuss 
how this mechanism could be realized in the highly dynamic environment of 
growing tissues in multi-cellular organisms. 

Key words: Morphogenesis; Pattern formation; Gene regulatory networks; 
Positional information; Proportion regulation 



1. Introduction 



Understanding the molecular machinery that regulates development of mul- 
ticellular organisms is among the most fascinating problems of modern science. 
Today, a growing experimental record about the regulatory mechanisms involved 
in development is accumulating , in particular i n wel l- studied mode l-organisms 



as. 



e.g., Drosophila or Hydra (iTechnau et all . l2000l : iBoschl . |2003| ). StiU, the 



genomic details known today are not sufficient to derive dynamical models of 
developmental gene regulation processes in full detail. Phenomenological models 
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of developmental processes, on the other hand, are well established today. Pio- 
neering work i n this field was done by Turing, who in his seminal paper in 1952 
(| Turing . 1X9521 ) considered a purely physico-chemical origin of biological pattern 
formation. His theory is based on an instability in a system of coupled reaction- 
diffusion equations. In this type of model, for certain parameter choices, stochas- 
tic fluctuations in the initial conditions can lead to self-organization and mainte- 
nance of spatial patterns, e.g. concentration gradients or periodic patterns. This 
principle has been successfully applied to biological morphogenesi s in nu merous 
applications ( Gierer and Meinhardtj . 19721 : Meinhardt and Gierei . 2000l ). How- 
ever, as experiments make us wonder about the astonishingly high com plexity 



20021) . they 



of single regulating genes in development (jBosch and Khalturinl . 
also seem to suggest that diffusion models will not be able to capture all de- 
tails of developmental regulation, and point at a complex network of regulating 
interactions instead. 

In theoretical work on pattern forma tion, both t he crucial role of local, 
induction-like phenomena in development (jSlackl . 1993) and limitations of diffu- 
sion based mechanisms in cellular environments ( Reillv and Melton , 19961 ) have 
lead to consideration of models that rely on local signal transfer via membrane- 
bound receptor s . A well-studied model in this context is juxtacrine signaling 
( Wearing et al. . 2000l : Owen et all . l2000l ). Positive feed-back, combined with 
juxtacrine signaling, can lead to gene ration of spatial pa. tterns with wavelengths 



that extend over many ceU lengths (jOwen et all . l2000t) . Further, it has been 



shown that a relay mechanism based on juxtacrine signaling can also lead to 
travelling wave fronts, a nd hence pro vide an alternative mechanism of long- 
range pattern regulation ( Monkl . 19981 ). The emergence of sharp spatial expres- 
sion boundaries of many genes in development, besides other genes that exhibit 
more graded profiles, is notoriously hard to explain in reaction-diffusion-based 
models. Recently it was shown in models of homeoprotein intercellular transfer 
( Kasatkin et al.l . [2007t iHolcman et al.l . 2007 ) that restricted local diffusion of a 
morphogen regulating its own expression can generate a morphogenetic gradi- 
ent. When two of t hese gradients meet, f or certain parameter values a sharp 



boundary is created (jKasatkin et al.l . 120071 ) 



The role of information processing in gene regulatory networks during devel- 
opment has entered the foc us of theoretical rese arch only recently. One pioneer- 
ing study was published by Jackson et al.l ( 198d ). who investigated the dynamics 
of spatial pattern formation in a system of locally coupled, identical dynamical 
networks. In this model, gene regulatory dynamics is approximated by Boolean 
networks with a subset of nodes communicating not only with nodes in the (in- 
tracellular) network, but also with some nodes in the neighboring cells. Boolean 
networks are minimal models of information processing in network structures 
and have been di s cussed as models of gene regulation since the end of the 1960s 
(jKauffman . Il969l 119931 : iGlasd . Il973l) . The model of Jackso n et al. demonstrated 
the enormous pattern forming potential of local information processing. More 
elaborate models that include a Boolean network description of cell-internal 
gene regulatory networ ks, local induct ive inter-cellular signals and a discrete 
model of cell adhesion (jHogeweg . l200nf ) point at a complex interplay between 
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regulatory dynamics, cell differentiation and morphogenesis. 

Along similar lines, Salazar- Ciudad et al. introduced a gene network mode l 
based on continuous dynamics (jSalazar-Ciudad et al. , I2OOOI : ISole et all . l2002h 
and coupling their networks by direct contact induction. Interestingly, they 
observe a larger variety of spatial patterns than Turing-type models with diffu- 
sive morphogens, and find that patterns are less sensitive to initial conditions, 
with more time-independent (stationary) patterns. This matches well the intu- 
ition that networks of regulators have the potential for more general dynamical 
mechanisms than diffusion driven models. Furthermore, dynamical models of 
regulatory networks that control basic st ages in development, e.g . the segment 
polarity network in Drosophila embryos ( von Dassow et all 120001) . have shown 
that d evelopmental module s are e xtremely robust against large parameter vari- 
ations; |AlberF^^^^Othmer ( 20031 ) even showed that the topology of regulatory 
interactions alone in a Boolean network model is sufficient to correctly predict 
both wild type and mutant patterns generated by the segment polarity net- 
work. Considering the temporal succession of regulatory dynamics rather than 
spatial patterns, similar results were obtained for other gene regulatory net- 
works import ant for cell development, e.g. the cell cycle n etwork of different 
yeast species (|Li et al.l . |2004 iDavidich and Bornholdtl . [2007h . The fact that, in 
many instances, simple discrete dynamical network models are sufficient to cap- 
ture essential properties of developmental dynamics, suggests that information- 
transfer-based processes controlled by the topology of regulatory interactions 
both within and between cells are very important for the extreme robustness 
and reliability observed in development despite considerable noise and large re- 
arrangements of cell ensembles due to cell proliferation and -movements. In this 
paper, we follow this new paradigm of interacting networks in pattern formation 
and in particular consider information-transfer-based processes. 

We start with a particular problem of position-dependent gene activation, 
as motivated from similar observations in Hydra. This animal is one of the most 
basal metazoa and exhibits extremely precise regulation of expression bound- 
aries under continuous cell movements. Also, it has remarkable properties to 
regenerate de novo after dissociation of cells, and to regulate its body pro- 
portions during growth. We introduce a novel, two-level theoretical approach 
to model pattern formation problems of this type: First, a coarse-grained de- 
scription in a deterministic cellular automata model is developed, which then 
is extended to a detailed model based on locally coupled, discrete dynamical 
networks. We show that this deterministic model explains both de-novo pat- 
tern formation after randomization of the pattern, and proportion regulation of 
gene activity domains. A threshold network model is derived which yields an 
upper estimate of the complexity of the regulatory module needed to solve the 
pattern formation task. Next, model dynamics is studied under noise and cell 
movements, and solved analytically in a mean-field approximation. It is shown 
that local, stochastic changes in gene expression states do not disrupt the spa- 
tial pattern, but contribute to its stabilization in the presence of cell flow by 
production of traveling domain boundaries (" quasi-particles" ) that coordinate 
global positional information, bearing some similarity to traveling waves found 
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in models based on juxtacrine signaling (jMonld . Il998f ). This suggests an inter- 
esting, new mechanism for reliable morphogenetic control, that might well apply 
to different types of tissues with high demands on regeneration. It is shown that 
the mechanism also works in growing tissues, explains pattern restoration after 
cutting the tissue in half, and is robust against noise in the detection of the 
body axis direction. Last, potential applications of the model are discussed. 



2. Motivation and model 

Let us first introduce and define the morphogenetic problem that we will use 
as a motivation for our novel pattern formation model. 

2.1. Proportion regulation in Hydra 

A classical model organism for studies of position dependent gene activation 
is the fresh water polyp Hydra, which has three distinct body regions - a head 
with mouth and tentacles, a body column and a foot region. The positions of 
these regions are accurately regulated along the body axis (Fig.[T]). 

Hydra also has the capacity to reproduce asexually by exporting surplus 
cells into buds; again, the position along the body axis where budding occurs 
is pr ecisely regulated at a bout two-thirds of the distance from the head to the 
foot ( Schiliro et al. . 19991 ). A number of genes are involved in regulation of the 



foot region and the budding zone; for example, both Pedibin and CN-NK2 are 
only ex pressed in the foot reg ion and turned off approximately in the budding 



region (jThomsen et al.l . 120041 ). While the CN- NK2 express i on do main exhibits 



a rather graded decay in the budding region ( Grens et al.L Il996h . it has been 



observed that, for example. Hedgehog ( Hh) is turned off precisely ]nsi below the 
budding region with a sharp boundary ( Kaloulisl . 1200 01. The relative position 



of the budding region and of the gene expression domains below it is almost 
independent of the animal's size, i.e. the ratio «/(! — a) (as denoted in Fig. 
[1]) is almost invariant under changes of body size. As the Pedibin/ CN-NK2 
system presumably plays an important role in determining the foot region, this 
invariance appears to be an essential prerequisite for maintaining the correct 
body proportions (proportion regulation) and to establish the head-foot polarity. 
Very precise regulation over a 10-fold size range has also been reported for the 
head- body proportion in Hydra, with a value close to 1/3 ( Bode and Bodd . 



Il980f l. Such precise regulation of position informat ion and body p roportions is 
a quite general problem in biological development ( Wolpertl . Il969l ). 



An interesting problem is how the specific properties of this regulation can 
be achieved by a small network of regulatory genes and if so, whether local 
communication between the cells (networks) is sufficient. This basic question 
is the central motivation for the present study. In particular, we consider the 
simplified problem of regulating one domain as, for example, the foot region 
versus the rest of the body. We consider this as a one dimensional problem as 
first approximation to the well-defined head-foot-axis in Hydra. We should note, 
however, that we do not intend to model in detail the regulatory mechanisms 
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underlying Hydra pattern formation, and rather take the observations made for 
this basic metazoon as an inspiration to introduce and study a simple, generic 
model of pattern formation. 



head 



1-a 



body 



foot 




Figure 1: Gene expression domains in Hydra, here for the example of the "foot" genes 
Pedibin (Ped) and Cn-NK2, and Hedgehog (Hh). Ped and Hh expression are bounded 
towards the body region of the animal; while Ped exhibits a graded decay in the 
budding region, Hh exhibits a sharp boundary. The relative position of the budding 
region and the associated expression boundaries, given by the ratio q/(1 — a), is 
independent of the absolute size of the animal (proportion regulation). Details are 
explained in the text. 



Developmental processes exhibit an astonishing robustness. This often in- 
cludes the ability of de novo pattern formation, e.g., to regenerate a Hydra even 



after complete dissociation of the cell ensemble in a centrifuge (jGierer et al 
Further, they are robust in the face of a steady cell flux: Hydra cells 



constantly move from the central body region along the body axis towards the 
top and bottom, where they differentiate into the respective cell types according 
to their position on the head-foot axis. The global pattern of gene activity is 
maintained in this dynamic environment. Let us take these observations as a 
starting point for a detailed study how the interplay of noise- induced regulatory 
dynamics and cell flow may stabilize a developmental system. 
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Figure 2: Diagram showing the interaction structure of the minimal network needed to 
solve the asymmetric expression task. For the sake of clarity, intracellular interactions 
between the two genes Gi and G2 are shown only for cell i, and likewise outgoing 
intercellular signals from the two genes two the neighbor cells i — 1 and i + 1 were left 
out. The transcription factors produced by gene Gi and G2 in cell i — 1 couple to the 
receptor systems Ri and R2, respectively, whereas in cell i + 1 the transcription factors 
produced by these genes couple to the receptor systems R3 and R4 (biased signaling). 
In cell i, the receptors release factors which regulate the activity of Gi and G2. 



2.2. One dimensional cellular automata: Definitions 

We here undertake a three-step approach to find a genetic network model 
that solves the pattern formation problem outlined above. In the first step, 
wc sum marize the properti es of the cellular automata model introduced in 
((Rohlf and Bornholdtl . [2005h . Cellular automata as dynamical systems discrete 
in time and stat e space are known to display a wide variety of complex patterns 
(LWolfra4 Il98l Il984allbh and are capable of solving complex computational 
tasks, including universal computation. We searched for solutions (i.e. rule 
tables which solve the problem) by aid of a genetic algorithm (for details, see 
Appendix). Candidate solutions have to fulfill four demands: Their update dy- 
namics has to generate a spatial pattern which 1) obeys a predefined scaling 
ratio a /{I — a), 2) is independent of the initial condition chosen at random, 
3) is independent of the system size (i.e. the number of cells Nc) and 4) is 
stationary (a fixed point). In the second step, this cellular automata rule table 
is translated into (spatially coupled) Boolean networks, using binary coding of 
the cellular automata states. The logical structure of the obtained network is 
reduced to a minimal form, and then, in step three, translated into a threshold 
network. 

To define a model syste m that performs the pattern formation task of domain 
self-organization (Rohlf and Bornholdd . 20051). conside r a one-dimensional cel- 
lular automaton with parallel update (jWolframl . Il984al) . Nc cells are arranged 
on a one-dimensional lattice, and each cell is labeled uniquely with an index 
i G {0, 1, Nc — 1}. Each cell can take n possible states S {0, 1, .., n}. The 
state (7i{t) of cell j is a function of its own state (Ti{t — 1) and of its neighbor's 



7 



states (Ti^i{t — 1) and ai+i{t — 1) at time t — 1, i.e. 



(7,(t) = / - 1), (T,(t - 1), (7,+ l(i - 1)] 



(1) 



with / : {0, 1, n}^ t—i- {0, 1, n} (a cellular automaton with neighborhood 3). 
At the system boundaries, we set cr_i = ajsi^^. = const. = 0. Other choices, e.g. 
asymmetric boundaries with cell update depending only on the inner neighbor 
cell, lead to similar results. The state evolution of course strongly depends 
on the choice of /: for a three-state cellular automaton {n — 3), there are 
3^^ ~ 7.626 • 10^^ possible update rules, each of which has a unique set of 
dynamical attractors. As we will show in the results section, n = 3 is the 
minimal number of states necessary to solve the pattern formation problem 
formulated above. 

Now we can formulate the problem we intend to solve as follows; Find a set 
!F of functions (update rules) which, given an initial vector a — (ctq, cat^-i) 
sampled randomly from the set of all possible state vectors, within T update 
steps evolves the system's dynamics to a fixed point attractor with the property: 



where [x] yields the largest integer value smaller or equal to the argument x (this 
is needed since the product a ■ Nc may lead to non- integer values). The scaling 
parameter a may take any value < a < 1. For simplicity, it is fixed here to 
a = 0.3. Notice that a does not depend on Nc, i.e. we are looking for a set of 
solutions where the ratio of the domain sizes r :— a/{l — a) is constant under 
changes of the system size, as motivated in section 12.11 by similar observations 
of proportion regulation in developing multi-cellular organisms. This clearly is 
a non-trivial task when only local information transfer is allowed. The ratio r 
is a global property of the system, which has to emerge from purely local (next 
neighbor) interactions between the cell's states. 

2. 3. Pattern formation by locally coupled Boolean networks 

One can now take a step further towards biological systems, by transferring 
the dynamics we found for a cellular automata chain onto cells in a line that 
communicate with each other, similar to biological cells. Ident ifying different 
dynamical states with (differentiated) cell types (Kauffmaj, Il969l ) and assuming 
that all model cells have an identical network of regulators inside, each of them 
capable to reproduce the rules of a cellular automaton by means of a dynamical 
coupling between subsets of regulators in direct neighbor cells, we obtain a model 
mimicking basic properties of a biological genetic network in development. 

Cellular automata rule tables can easily be formalized as logical (Boolean) 
update tables, e.g., for n = 3, two internal nodes with states (j\,cr\ G {0, 1} can 
be used for binary coding of the cell states (where i labels the cell position). 
One then has 





{<jl^{t),am} - fi,2<yl:2\t-l),al^{t - l),al+;{t-l)) 



(3) 
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with /i,2 : {0, 1}^ {0, 1}^. The so obtained rule tables are, by application 
of Boolean logic, transformed into a minimized conjunctive normal form, which 
only makes use of the the three logical operators NOT, AND and OR with a 
minimal number of AND operations. This is a rather realistic approximation for 
gene regulatory networks, as the AND operation is more difficult to realize on the 
basis of interactions between transcription factors. Other logical fu nctions as, 
e.g., XOR, are even harder to realize biochemically (jDavidsonl . 120011 ). However, 
we should notice again that it is not our intention to develop a specific or 
detailed biochemical model for the proposed pattern formation problem, but 
rather to find the simplest possible network model that reproduces the basic 
phenomenology. Nonetheless, a model of this type may well serve as a starting 
point for subsequent, more sophisticated models. The basic structure of the 
constructed network and a possible biological interpretation is shown in Fig. [2l 
Concordant with the spirit we followed so far, we assume that symmetry of signal 
transmission is broken on a local-cell scale, without specifying this in detail; 
several self-organizi ng mechanisms are conceivabl e in cellular systems, e.g., local 
chemical gradients (jGurdon and Bourillot . 12001 ) or anisotropic distribution of 



receptors on cell membranes ( Galle et al. . 20021 ). In Fig. 2, for the sake of 



simplicity, the specific example of a biased distribution of signal-transmitting 
receptors was chosen. 



i- 1 



input layer M 



hidden layer 



pattern gene 




i + 1 



a[=sign(Ec|,.ai +h) 



Figure 3: Schematic sketch of a threshold network controlling spatial gene activity 
patterns. Signals from genes in direct neighbor cells constitute the input layer of the 
network. In the hidden layer, this information is processed, with logical functions 
implemented as weighted sums of the inputs. The output of the genes in this layer 
then controls the pattern genes. Notice that there is feed-back to the hidden layer, as 
well as to the neighbor cells (dashed arrows). 



2.4- The simplest dynamics: locally coupled threshold networks 

Perhaps the simplest dynamical model for transcriptional regulation net- 
works are threshold networks, a subset of Boolean networks, where logical 
functions are modeled by weighted sums of the nodes' input states plus a 
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threshold h (jKurtenl . Il988f ). They have proven to be valuable tools to ad- 
dress questions associate d to the dynamics and e volut i on of gene regulatory 
netw orks (Wagner . 1994: Bornholdt and Sneppen , I2OOOI: iBornholdt and Rohli 
l200d iRohlf and Bornholdtl l2002l hOoH) . 

Any Boolean network can be coded as a dynamical threshold network with 
suitable thresholds assigned to each network node. For the system of coupled 
networks discussed in this paper, this network contains minimally three hier- 
archies of information processing ( "input layer" : signals from the genes in the 
neighbor cells at time t ~ 1, "hidden layer": logical processing of these signals, 
"output layer" : states of the two "pattern genes" al and ah in cell i at time t 
(see Fig. [3] for a schematic sketch of the system structure) The genes' states 
now may take values ct^ ~ ±1, and likewise for the interaction weights one has 



±1 for activating and inhibiting regulation, respectively, and 



if 



gene i does not receive an input from gene j in cell I. The dynamics then is 
defined as 

a;.(t)-sign(/,(t-l)) (4) 

with 



(5) 



fc=l 



for the "hidden" genes (compare Fig. [S]) , where al,k S {1,2} is the state of 
the fcth pattern gene in cell I (there are no couplings between the genes in the 
"hidden" layer). The threshold hj is given by 



i+l 



-k] I 



(6) 



k=l l=i-l 



which implements a logical OR operation. For the "output" (pattern) genes Gi 
and G2 in cell i, one simply has 



fkit) 



E 

1=1 



(7) 



i.e. the weights are all set to one, and the (negative) threshold equals the number 
of inputs kf^ that gene Gk,k G {1,2} receives from the hidden layer genes 
(logical AND). 

3. Results for deterministic dynamics 



Let us briefly summarize the dynamics of the simple stochastic cellular au- 
tomata model of spatial pattern formation based on local information transfer 



^Notice that this structure is quite similar to a feed-forward neural network, however, in 
our system there is regulatory feed-back from the output-layer to the 'hidden' layer. 
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( Rohlf and Bornholdtl . and its de novo pattern formation by generating 

and regulating a domain boundary. Subsequently, we will discuss in detail how 
this very general mechanism can emerge as a result of interacting nodes in cou- 
pled identical networks, as a model for gene regulation networks in interacting 
cells. 



3.1. Cellular Automata Model 

The first major outcome of the cellular automata model is that a number 
of n = 3 different states is necessary and sufficient for this class of systems 
to solve the given pattern formation task 0. The update table of the fittest 
solution found during optimization runs, which solves the problem independent 
of system size for about 98 percent of (randomly chosen) initial conditions (i.e. 
has fitness $ = 0.98), is shown in Table I. 




Figure 4: A typical dynamical run for the automata as defined in Table I, here for 
a system of size Nc ~ 250 cells (deterministic dynamics, no noise), starting from a 
random initial configuration. Time is running on the y-axis from top to bottom. Cells 
with state ai = are depicted in black color, cells with tJi = 1 in red and cells with 
Ui = 2 in blue. 



^The case n = 2 corresponds to the class of elementary (Wolfram) cellular automata with 
a very restricted set of 256 possible update rules. In our extensive genetic algorithm runs, no 
solution for the here considered problem was found for n = 2, even under diverse variations 
of the boundary conditions 
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Figure 5: Finite size scaling of the self-organized relative domain size a as a function of 
the total number of cells Nc- In the limit of large system sizes, a converges towards a 
fixed value Ooo = 0.281 ± 0.001 (as denoted by the straight line fit). The inset shows 
the finite size scaling of the variance Var[a{Nc)]', the straight line in this log-log plot 
has slope —1.3 and indicates that fluctuations vanish with a power of the system size. 

Fig. m shows the typical update dynamics of this solution. The finite size 
scahng of the self-organized relative domain size a as a function of the number 
of cells Nc is shown in Fig. [S] In the limit of large system sizes, a converges 
towards 

ttoo 0.281 ±0.001. (8) 

The variance of a vanishes with a power of Nc, i.e. the relative size of fluc- 
tuations induced by different initial conditions becomes arbitrarily small with 
increasing system size. Hence, the pattern self-organization in this system ex- 
hibits considerable robustness against fluctuations in the initial conditions. The 
main mechanism leading to stabilization at Ofoo = 0.281 is a modulation of 
the traveling velocity Vr of the right phase boundary in Fig. [4] such that the 
boundary on average moves slightly less than one cell to the left per update 
step, whereas the left boundary moves one cell to the right exactly every third 
update step {vi = 1/3). The modulation of the right boundary can be seen as 
the result of interacting phase boundaries reminiscent of particle interactions. 
This pict ure of "particle computat i on" is a useful concept also in various other 
contexts (jCrutchfleld and Mitchell . Il995l ). 



From the fact that cells interact only with nearest neighbors one might con- 
clude that three cells in a row in principle would be sufficient to generate a 
pattern, which would be in clear contradiction with a substantially larger min- 
imum size of aggrega tes that was found, for example, in the case of Hydra 
( Technau et al. . 2000) ■ The results summarized in Fig. [5l however, indicate 
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Table 1: Rule table of most successful cellular automata solution found during genetic 
algorithm search. In the left column, the rule table index is shown, running from to 
26, in the middle column the three input states at time t are shown, the right column 
shows the corresponding output states at time t 

that pattern formation becomes very inprecise for systems smaller than 100 
cells and typically fails for N < 25. Hence, the proposed mechanism is compat- 
ible with a required minimum system size that substantially exceeds the range 
of local communication. While in reaction-diffusion based models of pattern 
formation a certain extension of the field is required in order that the different 
diffusion rates come into play, in our model the differential propagation of phase 
boundaries leads to a similar effect. 

Let us now derive a quantitative model that approximates the system dy- 
namics. Since the left phase boundary travels at a constant speed of vi = 1/3, 
we only have to derive a (stochastic) model for the absolute value | (vr) \ of the 



configuration 


PijO 




(2,1,0) 


0.1646 


0.257 


(0,2,0) 


0.1852 


0.591 


(2,2,0) 


0.2058 


0.167 


(1,2,0) 


0.1646 


0.862 


(1,1,0) 


0.1317 


0.778 


(0,1,0) 


0.1482 


2.629 



Table 2: The six possible configurations at the right phase boundary with their respec- 
tive probabilities Pijo and velocity contributions (u)ijo of the corresponding transition 
trees (compare Fig. [6]). 
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X,(2,l,0) 




Figure 6: Transition tree for the boundary configuration (2, 1,0). Depending on the 
state of the left cell X, transition to different configurations occur. Numbers on arrows 
indicate the total number of the respective branches, the numbers at the bottom are 
the velocity (boundary readjustment) contributions of the respective branches. For 
details, see text. 



average traveling speed of the right phase boundary; the equilibrium boundary 
position then follows as 

"=^^T • (9) 

\{Vr)\+Vl ^ ' 

At the right phase boundary, there are three configurations (local update neigh- 
borhoods) that do not lead to a readjustment of the boundary (namely, (0, 2, 0), 
(2, 1,0) and (2,2,0), the zero on the right marks the boundary). The configu- 
rations (1,2,0) and (1, 1,0) readjust the boundary one cell to the left, whereas 
the (extended) configurations (xi, 0:2, 0, 1, 0) move the boundary either two or 
three cells leftward, depending on the states x\ and Xi. Using a Markovian ap- 
proximation (i.e., a one-step master equation neglecting transition correlations 
between the six boundary configurations), \{Vj)\ is approximated by 

KMI =0-Po + l-Pi + 2-p2 + 3-p3, (10) 

where pj are the respective probabilities to have a configuration that leads to 
boundary readjustment i cells at the left at the next time step. We neglect the 
slight asymmetries in the rule table and assume that each state a G {0,1,2} 
appears with probability 1/3, hence it is straight-forward to derive po = 1/2 
and pi = 1/3. A slightly more detailed analysis yields p2 = 1/18 and ps = 1/9, 
leading to 

Kv^)i| «0.78, (11) 
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which is about 9% below the true value \{vr)\ ~ 0.852 measured in model 
simulations. 

To improve the approximation, we now take into account transition correla- 
tions between different configurations and the slightly asymmetric state distri- 
bution in the rule table. For each of the six boundary configurations, a transition 
tree similar to Fig. [S]is derived (for the last configuration, this tree consists of 
only one time step and one transition, i.e. collapses on the Markovian approxi- 
mation). Taking the average over the velocity contributions V2i of all branches 
of the second time step (where i numbers the branches), the contribution of the 
whole tree to the average phase boundary velocity per time step is 

1 

{v)tree = 7^ Y] W2j, (12) 

where is the number of branches (here, n2 — 9). The start configurations 
with their respective probabilities p^o and velocity contributions {v)ijo of the 
corresponding transition trees are listed in table 1. The phase boundary velocity 
now is calculated as the weighted average 

2 

\{Vr)\= P^JO ■\{v)tjo\- (13) 

Inserting the values of table [2l one finds 

\{Vr)2\~0-S2. (14) 

Obviously, this value is a much better estimate than the zero-order approxima- 
tion {vr)i, but still 4% below the value \ {vr)\ « 0.852 measured in simulations. 
We conclude that this difference is an effect of higher order correlations not 
included in our analysis. 

3. 2. Interaction topology of the minimal network 

In this section, let us derive the structure of a minimal Boolean network that 
solves the pattern formation problem, based on the previously discussed results 
for cellular automata. We will see that this network has biologically realistic 
properties regarding the number of genes necessary for information processing 
and the complexity of interaction structure, making it well conceivable that 
similar "developmental modules" exist in biological systems. 

3.2.1. Boolean network model 

Let us now undertake the first step from the previous, coarse-grained model 
of pattern formation to a detailed model that takes into account the information 
processing capacity of cell- internal regulatory networks, that can communicate 
locally with neighboring cells. The rule table summarized in Table 1 is easily 
formalized in binary coding, i.e 00, 1 ^ 01 and 2 10, this corresponds 
to two "genes" Gi and G2 one of which (Gi) is active only in a domain at the 
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Figure 7: Boolean representation of the minimal network, minimized conjunctive nor- 
mal form. G\ with a G {1,2} and b £ {i — + 1} denotes gene a in cell number 
b. The inputs in the left branches of the trees are given by the genes' states at time 
t - 1. ^ denotes NOT, denotes logical AND and ® logical OR. 



left side of the cell chain. The so obtained Boolean update table is reduced 
to its minim i zed co njunctive normal form, using a Quine-McCluskey algorithm 
l McCluskevLll956h . ?or the construction of the network topology we use the 
conjunctive normal form, as it is a somewhat biologically plausible solution 
with a minimal number of logical AND operations. In principle, other network 
topologies, e.g. with more levels of hierarchy, are possible and biologically plau- 
sible, however, they involve a higher number of logical sub-processing steps, i.e. 
a higher number of genes, hence we will not discuss them here. 

Considering the huge number of possible input configurations which the 
outputs theoretically could depend on, the complexity of the resulting network 
is surprisingly low. As shown in Fig. [71 the output state of gene Gi only depends 
on five different input configurations of at maximum four different inputs, gene 
number two on six different input configurations of at maximum four different 
inputs. This indicates that the spatial information flowing into that network 
is strongly reduced by internal information processing (only a small number of 
input states leads to output "1"), as expected for the simple stationary target 
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Figure 8: Threshold network realization of the pattern formation system. Solid line 
arrows denote links with Wij = +1, dashed arrows denote hnks with Wij = — 1. The 
inputs from the genes in the neighbor cells {G]~^ ,G'2~^ , G^^^^ and G^^^^) are processed 
by a layer of "hidden genes" (colored circles in the middle of the scheme) with different 
thresholds h implementing a logical OR operation on the inputs. The processed signals 
then are propagated to the two pattern genes Gi and G2. The threshold of gene Gi 
is hi — —5, for gene G2 one has h2 = —6 (logical AND). Notice several feed-back 
connections from genes G\ and G2 to the hidden layer. 

pattern. Nevertheless, this information processing is sufficient to solve the non- 
trivial task of domain size scahng. 

3.2.2. Coupled threshold network model 

Threshold dependence of the states of regulatory elements constitutes a bio- 
chemically simpler paradigm of switching behavior; information processing dy- 
namics is encoded in activating and inhibiting interactions only, without the 
need for complex Boolean update tables. The simpler switching dynamics comes 
at the expense of an increased network size, hence the formalization as a thresh- 
old network gives us an estimate for the upper limit of regulatory network size 
needed to solve the pattern formation problem. The coupled threshold network 
system, that was derived according to the method outlined in section 12. 4[ is 
shown in Fig. [8] The states of the genes Gi and G2 at time t in a cell i and its 
two neighbor cells i — 1 and i -I- 1 serve as inputs of 11 information-processing 
genes ("hidden" layer). The state of these genes then defines the state of Gi 
and G2 in cell i at time t + 2 (output layer) . Additionally, there is some feed- 
back from Gi and G2 to the information processing layer, as expected for the 
dependence on cell-internal dynamics already present in the cellular automata 
implementation of the model. 

The resulting stationary spatial patterns of the information-processing 'hid- 
den' genes and of genes Gi (the "domain gene") and G2 (active only at the 
domain boundary) are shown in Fig. [5] (snapshots of five different update time 
steps for a system of 70 cells). Starting from a random initialization of the 
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t = 40 




Figure 9: Snapshots of spatial gene activity patterns of the network shown in Fig. 6., 
for a system of 80 coupled networks (cells) at different update times t. Each row shows 
the state of one network gene, color coding is the same as in Fig. |8l if the respective 
gene is not active in the cell at time t, the cell is shown in black color. The two bottom 
rows show the states of the two target (pattern) genes. After 55 update time steps, 
the target pattern (compare Fig. [4]) has self-organized. 



two pattern genes Gi and G2, due to the high-level genetic information pro- 
cessing in the hidden layer the target pattern self-organizes robustly within 55 
update time steps. The network we construct here, regarded as a "develop- 
mental module" defining the head-foot polarity through spatially asymmetric 
gene expression, has a size similar to comparable biological modules (compare . 



for example, the segment polar ity network in Drosophila (jvon Dassow et al 



200d: [Albert and Othmerl . |200i)) as well as similar complexity (e.g., average 



connectivity K m 3). 



4. Dynamics under noise and cell flow 

In the following we will study dynamics and robustness of the model with re- 
spect to noise. Two kinds of perturbations frequently occur: Stochastic update 
errors and external forces induced by a directed cell flow due to cell prolifera- 
tions. Both types of perturbations are very common during animal development, 
e.g., in Hydra cells continuously move from the central body region along the 
body axis towards the top and bottom, and differentiate into the respective cell 
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Figure 10: Quasi-particles, started by stochastic update errors, lead to control of the 
boundary position under noise (left panel). The F particle (top right) leads to read- 
justment of the boundary two cells to the left, the A particle (bottom right) leads to 
readjustment of the boundary one cell to the right. 

types along the way according to their position on the head-foot axis. However, 
similar problems of reliable pattern formation in noisy and highly dynamical 
environments occur in tissues with a high turnover. High proliferation rates 
and/or movement of cells occur, for example, in skin tissues and when stochas- 
tic phenomena of cell (re-)differentiation, e.g. stochastic stem cell production, 
are found. In our model, we abstract stochastic changes in differentiation states 
by stochastic update errors, and consider directed cell movements in the form 
of a steady cell flow. 

Let us define stochastic update errors with probability p per cell, leading 
to an average error rate re — pNc- Interestingly, this stochastic noise starts 
moving "particle" excitations in the cellular automaton which, as a result indeed 



Figure 11: For moderate error rates re, the domain boundary is stabilized at an average 
position a* — 1/3 (left panel, — 0.1). Around « 0.2, there is a crossover to a 
domain size vanishing with (middle panel). In the right panel, the high noise limit 
is shown, with a considerably shrinked blue domain due to strong particle interference 
{re = 2.0). 
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Figure 12: Average boundary position Q* as a function of the error rate for system 
sizes Nc = 100, Nc = 400, and Nc = 1600. Tfie abscissa is logaritfimic. Numerical 
data are averaged over 200 different initial conditions with 2 • 10^ updates each. The 
dashed curves show the mean field approximation given by Eqn. (|35() . the straight 
dashed lines mark the unperturbed solution a* = 0.281 and the solution under noise, 
a* = 1/3. 



stabilize the developmental structure of the system. To prepare for the details 
of these effects, define first how we measure the boundary position properly in 
the presence of noise. Let us use a statistical method to measure the boundary 
position in order to get conclusive results also for high p: Starting at i = 0, we 
put a "measuring frame" of size w over cell i and the next w — 1 cells, move 
this frame to the right and, for each i, measure the fraction z of cells with state 
(7 = 2 within the frame. The algorithm stops when z drops below 1/2 and the 
boundary position is defined to be i + 'w/2. 

One can show that, for not too high p, there are only two different quasi- 
particles (i.e. state perturbations moving through the homogeneous phases), as 
shown in Fig. [TOl In the following, these particles are called T and A. The 
r particle is started in the (72 phase by a stochastic error ai = 2 ai ^ 2 at 
some i < aNc), moves to the right and, when reaching the domain boundary, 
readjusts it two cells to the left of its original position. The A particle is started 
in the cto phase by a stochastic error = — s- (Ti 7^ at some i > aNc and 
moves to the left. Interaction with the domain boundary readjusts it one cell to 
the right. Thus we find that the average position a* of the boundary is given 
by the rate equation 

2Q*re = (l-a*)re, (15) 

i.e. a* = 1/3. Interestingly, for not too high error rates re, a* is independent 
from and thus from p. If we consider the average boundary position a* as a 
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Figure 13: For the system with stochastic update errors, fluctuations of the boundary 
position a around the average position q* = 1/3 are Gaussian distributed. The 
figure compares the numerically obtained stationary probability distribution with the 
analytic result of Eqn. (|23|l for three different system sizes. All data are gained for 
re = 0.1 and averaged over 100 different initial conditions with 2 ■ 10® updates each. 
The inset shows a typical timeseries of the boundary position. 



system-specific order parameter which is controlled by the two quasi-particles, 
then comparing the solution of Eqn. p5p to the equilibrium position in the 
noiseless case indicates that the system undergoes a step-like discontinuity with 
respect to a* ai p = = 0. This conclusion is supported by a numerical 
analysis of the finite size scaling of this transition (cf . appendix [C| . 

The solution a* = 1/3 is stable only for < rg < 1/2. As shown in Fig. [TUl 
the interaction of a T particle with the boundary needs only one update time 
step, whereas the boundary readjustment following a A particle interaction 
takes three update time steps. Hence, we conclude that the term on the right 
hand side of Eqn. |35)) . which gives the flow rate of A particles at the boundary, 
for large will saturate at 1/3, leading to 

2a*re = i (16) 

with the solution 

a*^\ + Q{Nc) (17) 

for > 1/2. Hence, there is a crossover from the solution a* = 1/3 to another 
solution vanishing with r~^ around rg = 1/2. The finite size scaling term Q(Nc) 
can be estimated from the following consideration: for p — > 1, the average 
domain size created by "pure chance" is given by a* = N(^^ '^n=o(^/'^)^ ' 
n « (3/4) A^p^. If the measuring window has size w, we obtain Q{Nc) ~ 



21 



0.01 



0.001 




— I r~ 

r, = 0.01 
analytic result 

r, = 0.005 
analytic result 

r, = 0.001 
analytic result 



0.0001 



200 



400 



600 



800 



1000 



Figure 14: Probability distribution p{t) of waiting times for boundary readjustments in 
the model with stochastic update errors for three different values of r^, semi-log plot. 
As expected for a Poisson process, p{t) is an exponential. 



{3/4:)wN^. To summarize, we find that the self-organized boundary position 
is given by 




0.281 ±0.001 if re = 

1/3 if 0<re<l/2 (18) 

(l/6)r-i + eiNc) if re > 1/2 



with a step-Hke discontinuity at re = and a crossover around = 1/2. 

Now let us consider the fluctuations of a around a* given by the master 
equation 

p^a) = 2arep'^-\a + 2d) + {l-a)rep'-\a-5) 

- {l~a)T,p-'-\a) (19) 

with 5 = 1/Nc- Eqn. ()19|) determines the probability p'^ (a) to find the boundary 
at position a at update time step r, given its position at time t — 1. This 
equation can be simplified as we are interested only in the stationary probability 
distribution of a. It is easy to see that the error rate r^ just provides a time 
scale for relaxation towards the stationary distribution and has no effect on 
the stationary distribution itself. Therefore, we may consider the limit — > 
^max .„ divide through r^. and neglect the last three terms on the right 
handside of Eqn. ([TO]) (which become zero in this limit). We obtain 

p^(a) = 2ap''-\a + 2d) + {I - a)p''-\a - d). (20) 
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Figure 15: Average domain size Q* as a function of the the cell flow rate r/ for three dif- 
ferent error rates Ve; numerical data (crosses and points) were sampled over 10 different 
initial conditions and le6 updates for each data point. "Left" indicates cell flow di- 
rected to the left system boundary, "right" to the right system boundary, respectively. 
The dashed lines are the corresponding solutions of Eqn. (|25p . 



To study this equation, we consider the continuum limit Nc — > oo. Let us 
introduce the scaling variables x = {a — a*)y/Nc, t = t/Nc and the probability 
density f{x,t) = Nc {a Nc)- Inserting these definitions into Eqn. ([20]) and 
ignoring all subdominant powers 0{1/Nc), we obtain a Fokker-Planck equation: 



dfix,t) 
dt 



dx^ 



3^a; ) f{x,t). 



The stationary solution of this equation is given by 



3 2 



fix) = \/ ^ exp 



(21) 



(22) 



i.e. in the long time limit t oo, the probability density for the boundary 
position a is a Gaussian with mean a* : 



p{a,Nc) 



3iVc 



27r 



■ exp 



iNc 



(a 



(23) 



From Eqn. we see that the variance of a vanishes ~ l/Nc and the relative 
boundary position becomes sharp in the limit of large system sizes. Fig. [T3] 
shows that this continuum approximation for Nc > 400 provides very good 
correspondence with the numerically obtained probability distributions. 

The stochastic nature of boundary stabilization under noise is also reflected 
by the probability distribution of waiting times t for boundary readjustments 
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due to particle interactions: the particle production is a Poisson process with 
the parameter X = r^, and the waiting time distribution is given by 

Pwait{t) = Te exp {-r^ t) (24) 

with an average waiting time {t) ~ r'^^. Fig. 1141 shows the waiting time distri- 
butions for different error rates re- 
in a biological organism, a pattern has to be robust not only with respect 
to dynamical noise, but also with respect, e.g., to "mechanical" perturbations. 
In Hydra, e.g., there is a steady flow of cells directed tow ards the animal's head 



and foot, due to continued proliferation of stem cells (jPavid and Campbelll . 



I972I ): the stationary pattern of gene activity is maintained is spite of this cell 



flow. Let us now study the robustness of the model with respect to this type 
of perturbation. Let us consider a constant cell flow with rate r/, which is 
directed towards the left or the right system boundary. In Eqn. (9), we now get 
an additional drift term on the left hand side: 



2aVe±r/ =re(l-a*), (25) 
with the solution 

[H'-'^) " (26) 

i if re < 

for the case of cell flow directed towards the left system boundary (plus sign 
in eqn. (|25|) ). One observes that a* undergoes a second order phase transition 
at the critical value r^'"'* = r/. Below r^"*^, the domain size a* vanishes, and 
above rg*"'* it grows until it reaches the value a'^ax = 1/3 of the system without 
cell flow. For cell flow directed towards the right system boundary (minus sign 
in eqn. ([25|l ). we obtain 

a* = | + 'f^'^^ (27) 

i 1 if r/>2re. 

In this case, the critical cell flow rate is given by rf — 2re, for cell flow rates 
larger than this value the CT2-domain extends over the whole system, i.e. a* = 1. 

Fig. [15] compares the results of numerical simulations with the mean field 
approximation of Eqn. (|25p . In numerical simulations, cell flow is realized by 
application of the translation operator Qai := Ui+i to all cells with < i < 
Nc — 1 every rj^ time steps and leaving umc-i unchanged. In case of cell flow 
directed to the right system boundary, in the limit r/ ^ 1 the boundary position 
a* detected in numerical simulations deviates from the mean field prediction, 
due to a boundary effect at the left system boundary (stochastic production of 
finite lifetime stationary oscillators, leading to intermittent flows of F particles 
through the system). 

To summarize this part, we see that in the model stochastic errors in dy- 
namical updates for r^ > rj indeed stabilize the global pattern against the 
mechanical stress of directed cell flow. 
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5. Proportion regulation in a growing system 



So far, we assumed that the system size N (the number of ceUs) is constant, 
which is a good approximation for an adult organism; in a developing organism, 
however, proportion regulation has to work under the condition of a steadily 
growing system size. Here, we study this problem for two simplified settings: 
first, for symmetric growth, i.e., new cells are added with probability 1/2 on 
either side of the chain of cells, and the growth rate Vg is constant on average; 



second, for homogeneous cell proliferation with probability j 
that daughter cells inherit the state of the mother cell. 



per cell, assuming 
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Figure 16: Asymptotic boundary position a^ir-g) in the case of unlimited symmetric 
growth at the boundaries of the cellular array, for four different error rates r^. 



5.1. Symmetric growth at the system boundaries 

Let us assume we start with a system of iVg cells, with an initial boundary 
position at cell iVi. In the deterministic case re = 0, it is straight-forward to 
see that the asymptotic boundary position in the limit of large times t is given 

by 

= lim a*{t) = ]- (28) 



(for details, see appendix ID. 1[) . This means that in the limit r^ = 0, proportion 
regulation cannot be maintained under the condition of a steady system growth. 
In the case Tg > and assuming infinite growth, the asymptotic boundary 
position is given by 

= lim a*{t) = lLl±I± (29) 

t-+oo To + 3re 
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Figure 17: Asymptotic domain boundary position aoo (lined curve) as a function of 
the initial system size A^'o, for homogeneous growth, as explained in the text. One 
has Qoo — ao — 1/No, where ao = 0.281 is the boundary position of the constant- 
size system (dashed line). Inset: Proliferation of the boundary cell (red) at time t — 1 
leads to readjustment of the boundary at its original position (indicated by the dashed 
line) at time t + 1, thereby increasing the black domain by one cell and hence slightly 
reducing a. 



(a derivation can be found in appendix ID.ip . Fig. 1161 shows a°°(rg) for four 
different values of r^', it becomes evident that an approximately 'correct' pro- 
portion regulation requires to be at the order of rg or larger, i.e. re/fg > 1- 
While Te (the rate of regulatory signals) may not be increased significantly above 
the growth rate rg, due to metabolic constraints, in later stages of development 
the steady decrease of rg will ensure that the condition r^/rg > 1 is fuUfilled 
and proportion regulation approaches the steady state of the adult organism. 

5.2. Homogeneous cell proliferation 

Another simple case is system growth by homogeneous cell proliferation. 
Assuming that daughter cells inherit the state of their mother cell, one can show 
that proportion regulation is maintained even in the case of zero noise, under 
the simplifying assumption that initial pattern formation takes place in a system 
of size A^o, and that system growth does not start before pattern formation has 
converged to its attractor. Due to an instability induced by proliferation events 
directly at the boundary cell with at = 1 (compare inset of Fig. \17\ . slight 
deviations asymptotic boundary position ao of non-growing systems are found 
for finite A'o (Fig. \17\ for details, cf. appendix ID. 2|) : 

ttoo = hm a{t) = ao - (30) 
t^oo 2Ao 

For the case when noise is present, it is not possible to find a general solution 
since proliferation can affect both the velocity and the type of particles travelling 
through the domains in intricate ways (in the case of symmetric growth at the 
system boundaries, as discussed in the previous subsection, this problem is 
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avoided). If pd is very small, however, we can assume that proliferation events 
and particle propagation are essentially decoupled and that the system has 
enough time to relax to a stationary state between proliferation events. In this 
limit, one can show that the asymptotic boundary position converges to the 
value a = 1/3 of the stationary size system as discussed in section 2] (for a 
derivation cf. appendix ID. 2p . 



6. Regeneration in a simulated cut experiment 

Simple multi-cellular organisms as, e.g., Hydra exhibit remarkable regenera- 
tion capacities, which include, as already discussed, proportion regulation and 
de-novo pattern formation after complete dissociation of the body tissue. Simi- 
larly, it was already observed in the late 19th century that polyps can be cut in 
half, leading to regeneration of two new, intact animals (). Without going into 
the more intricate details of these experiments, we now demonstrate that, given 
minimum level of noise in the system, our model in principle can reproduce this 
type of observation. Fig. (TS] illustrates a simulated cut experiment, where, after 
300 initial system updates, the cellular array was cut into two equal-sized halves. 
After just 500 subsequent updates, both new sub-systems have self-organized 
again into the target pattern with a = 1/3. 



7. Robustness under noisy direction recognition 

While asymmetries in receptor distribution on cell membranes, asymmet- 
ric distribution of cell factors in the cell or an extrenal gradient might provide 
some information about the asymmetry (the direction) of the spatial pattern 
along the body axis, which then can be processed by a cell- internal gene reg- 
ulatory network, a substantial amount of noise can be expected to be present 
in this process. In particular, in the system discussed in this paper, this type 
of information can be assessed only locally, hence we expect that local errors 
in 'direction recognition' can substantially disrupt the emergence of the global 
pattern. 

We now test the robustness of our model with respect to this type of errors. 
Let {ai^i{t) , (Ji{t) , ai+i{t)) be the state of the neighborhood of cell i at time t, 
then the state of cell i at time i -I- 1 is given by 



^ ./^^ ^-.^ f /(crt-i(i),cr,(i),cri+i(t)) with prob. 1 - Pdir /g^x 
\ f{(Ti+i{t),at{t),(7t-i{t)) with prob. pdir ' 

where pdir is the probability of false direction recognition, and /(.) is the cor- 
responding rule table entry associated to the state (cri-_i(t), (Ti(i), (Ti+i(i)) and 
its locally inverted state {ai+i{t), (7i{t), ai-i{t)), respectively. 

Our first finding is that the dynamics of the original system (deterministic 
dynamics, i.e. = 0), is indeed disrupted, due to a destabilization of the 
boundary state. For pdir < 0.2, a always goes to zero, while at pdir ~ 0.2, 
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Figure 18: A simulated cut experiment. After 300 updates, a system of 500 cells was 
cut into two sub-systems of 250 cells each (upper panel). After only 30 updates, 
in both sub-systems reorganization of the boundary position starts (middle panels). 
After about 500 system updates, both sub-systems have self-organized into the target 
pattern with q = 1/3 (bottom panels). In the simulation, = 0.01 was applied. 



there is an abrupt jump to a w 0.95 (Fig. [20l top left panel). However, the 
situation changes substantially in the much more realistic case of a finite error 
rate in dynamical updates, i.e. > 0. Fig. [12] demonstrates that in this 
case initial pattern formation (left panel), as well as control of the boundary 
position by noise induced particles (right panel) work, although the trajectory 
of the information-transmitting particles is blurred out and broadened by the 
stochastic errors in direction recognition. The latter effect is reflected by an 
increase in fluctuation size (an increase of the variance) of the boundary position 
a with increasing pdir, in particular in the limit of high dynamical error rate 
Te (Fig. [201 bottom panels). Remarkably, the average boundary position a is 
stabilized over a wide range of the new control parameter pdir , both in the limit 
of small Te (Fig. [20l left upper panel) and large Ve (Fig. [20l right upper panel). 
While there is some dependence on both pj^ij. and r^, as reflected by the fact 
that the curves a{pdir, fe) for different values of do not collapse, the principal 
pattern formation mechanism still works, given we stay at reasonable values 
Pdir < 0.2. Hence, the system exhibits remarkable robustness also with respect 



28 



Figure 19: Dynamics of pattern formation under noisy direction recognition; parameters 
in the simulation shown here were pdir = 0.1 and = 0.003. Initial pattern formation 
(left panel) as well as control of the boundary by quasi-particles (right panel) still 
work, though particle trajectories become broadened and blurred. 



to errors in direction recognition. Notice that this robustness was not selected 
for in GA runs, i.e., it is a truly emergent property of the system dynamics. In 
a real system of coupled gene regulatory networks, additional mechanisms for 
error correction might be present, that, for example, process information not 
only from direct neighbor cell s, or exploit the 2D or 3D geometry of a real tissue 
(|Rohlf and Bornholdtl . [2004bl ). 



8. Summary and Discussion 

In this paper we considered the dynamics of pattern formation motivated 
by animal morphogenesis and the largely observed participation of complex 
gene regulation networks in their coordination and control. We therefore chose 
a simple developmental problem to study toy models of interacting networks 
that control pattern formation and morphogenesis in a multicellular setting. 
In particular, the goal was to explore how networks can offer additional mech- 
anisms beyond the standard diffusion based process of the Turing instability. 
Our results suggest that main functions of morphogenesis can be performed by 
dynamical networks without relying on diffusive biochemical signals, but using 
local signaling between neighboring cells. This includes solving the problem of 
generating global position information from purely local interactions, but also 
it goes beyond diffusion based models as it offers solutions to developmental 
problems that are difficult for such models and avoids their inherent problem of 
fine tuned model parameters. Indeed, it has been shown in case studies that this 
paradigm applies well to development, as for example for the segment polarity 
network of Drosophila, which exhibits robustness again st parameter variations 
by several orders of magnitude ( von Dassow et al. , 200Clf ) , and where spatial gene 



expression patterns can be predicted reliably from the topology of regula tory 
interactions alone in a Boolean network model ([Albert and Othmer . |2003| ). In 



many cases, developmental processes as, e.g., the establishment of positional in- 
formation, may rely on this type of internal information processing rather than 
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Figure 20: Upper panels: average boundary position a as a function olpdir, for different 
values of r^. Left panel: low noise limit {r^ —> 0), right panel: high noise limit. Lower 
panels: Fluctuations (temporal variance) of the boundary position as a function of 

Pdir ■ 



on interpretation of global chemical gradients. In this type of local information 
processing, several ways hovif spatial symmetry of morphogenetic signals could 
be broken are concei vable. Cells potent ially could exploit local anisotropies in 
receptor localization teal le et al l. l2002h. as wel l as gr adients produced by lo- 
cal propagat ion of morphogens (Kasatkin et al. . 2007 ) or juxtacrine signaling 
( Monkl . Il998l ). In cither case, the mechanism proposed in our model would ex- 
hibit considerable robustness, since only a rough estimate on the direction of the 
receptor anistropy or the gradient is needed (compare section 7 on robustness 
under noisy direction recognition. 

The network model derived here performs accurate regulation of position in- 
formation and robust de novo pattern formation from random conditions, with a 
mechanism based on local information transfer rather than the Turing instabil- 
ity. Non-local information is transmitted through soliton-likc quasi-particles in- 
stead of long-range gradients. Two realizations as discrete dynamical networks, 
Boolean networks and threshold networks, have been developed. The resulting 
networks have size and complexity comparable to de velopmental gene regula- 



tion modules as observed in animals, e.g., Drosop hila (jvon Dassow et al. I. I2OOOI: 



Albert and Othmeil . [2003h or Hydra (|Boschl [2OO3I) . The threshold networks (as 



models for transcriptional regulation networks) process position information in a 
hierarchical manner; in the present study, hierarchy levels were limited to three, 
but realizations with more levels of hierarchy, i.e. more "pre-processing" of in- 
formation are also possible. Similar hierarchical and modul ar orga nization are 
typical signatures of gene regulatory networks in organisms (|Davidson. . ,200lh . 
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Robustness of the model was studied in detail for two types of perturba- 
tions, stochastic update errors (noise) and directed cell flow. A first order phase 
transition is observed for vanishing noise and a second order phase transition 
at increasing cell flow. Fluctuations of the noise-controlled boundary position 
were studied numerically for finite size systems and analytically in the contin- 
uum limit. We find that the relative size of fluctuations vanishes with 1/Nc, 
which means that the boundary position becomes sharp in the limit of large sys- 
tem sizes. This means that, based on the proposed local mechanism of coupled 
regulatory networks, positional information can be reliably controlled also in 
large tissues, which is problematic in the alternative case of morphogenetic gra- 
dients that are typically limited to relatively small spatial domains. Dynamics 
under cell flow was studied in detail numerically and analytically by a mean field 
approximation. A basic observation is that noise-induced perturbations act as 
quasi-particles that stabilize the pattern against the directed force of cell flow. 
Hence, we make the interesting observation that noise in local gene expression 
states (over several orders of magnitude in the relevant dynamical parameters) 
contributes to robustness of the global developmental dynamics; furthermore, 
this is a truly emergent property of the spatial system, which was not selected 
during simulated evolution. At a critical cell flow rate, there is a second order 
phase transition towards a vanishing domain size or a domain extending over 
the whole system, depending on the direction of cell flow, respectively. The 
proposed local mechanism of developmental pattern control also works in grow- 
ing tissues, reproduces pattern regeneration after cutting a tissue in half, and is 
robust against noise in the recognition of the body axis direction. 

Let us briefly compare the prospects and limitations our model with respect 
to other recent models suggested for pattern formation, and with experimen- 
tal evidence. "Local" models of pattern formation, in contrast to older models 
that require long-range diffusion (which is problematic in multi-cellular environ- 
ments in a nu mber of regards ), have been sugges ted in the context of juxtacrine 



signalling (JS. iMon k (19 9^); Owen et all ( 200C)) and homeopr otein intercellu- 

[itI 



lar transfer fHIT. [Kasatkin et all (200l1 ): lHolcman erall (|2007f )). Similar to JS 
models with relay, traveling waves/excitations emerge in our model as a means 
to provide long range communication. Sharpness and precision of boundary reg- 
ulation is shared with HIT models, where, however, this property arises from 
a different mechanism (meeting of morphogenetic gradients). While in HIT 
models noise can substantially affect boundary regulation, an essential prop- 
erty of the model analyzed in our study is its astonishing robustness against 
noise. When cell movement is present, noise in fact considerably contributes to 
pattern regulation and -stabilization. An evident limitation of the model arises 
from the fact that it accounts for regulation of sharp expression boundaries, but 
not for graded expression patterns. Sharp boundaries ar e indeed found for many 
genes in development (examples in Hydra are Hedgehog (iKaloulid . 20001 ) and the 

wh i le oth er genes 
20061) exhibit 



sharp basal bord er of HyBMP5-8b (Reinhardt et all. 120041)) 



such as CnNK2 ( Grens et al.l . 119961 ) and Dkk ( Augustin et al 



more graded expression patterns along the body axis. It seems quite natural to 
assume that, in addition to local mechanims as proposed in our model, other 
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mechanims of pattern formation are present in developing organisms that work 
on other scales and in different functional contexts, involving regulatory pro- 
cesses based on graded expression profiles. The hierarchical interplay of such 
diverse regulatory mechanisms might substantially contribute to the astonish- 
ing robustness of developmental processes. Going beyond basal metazoa such 
as Hydra, other interesting applications of our model are conceivable. For in- 
stance, local communication systems between adjacent cells as, for example, the 
Delta/Notch systems, play a decisive role in vertebrate development, with trav- 
eling waves providing long-ra nge synchronization of developmental processes 
(|Ozbudak and Pourgui^ . |2008| ). 

Several extensions of the model as described in this paper are conceivable. 
In the present model, the cell flow rate r/ is considered as a free parameter, the 
global pattern, however, can be controlled easily by an appropriate choice of the 
error rate re- This may suggest to extend the model by introduction of some 
kind of dynamical coupling between re and r/, treating r/ as a functioii of re . 
Interestingly, similar approa ches have been studied by Hogeweg ( Hogewea ^ 2000l) 
and Furusawa and Kaneko ( Furusawa . 2000t Furusawa and Kanekol 2003 ): In 
both models of morphogenesis, the rate of cell divisions is controlled by cell 
differentiation and cell-to-cell signaling. Dynamics in both models, however, is 
deterministic. An extension of our model as outlined above may open up for 
interesting studies how stochastic signaling events could control and stabilize a 
global expression pattern and cell flow as an integrated system. Other possible 
extensions of the model concern the dimensionality: In two or three dimensions 
other mechanisms of symmetry breaking might be present, possibly leading to 
new, interesting dynamical effects. 

A Java applet simulation of the model can be found at 



http : //www . theo-physik . uni-kiel . de/ ~rohlf /development . html 
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A. Genetic algorithm searches 

Let us briefly recapitulate here how the rule table of the model has been 
obtained by the aid of a genetic algorithm. 

A.l. Definition of the GA 

In order to find a set of update rules that solve the problem as formu- 
lated in section I I, cell ular automata have been evolved using genetic algorithms 
(|Mitchell et al.L Il994h . Genetic algorithms are population- based search algo- 
rithms, which are inspired by the interplay of random mutations and selection 
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as observed in biological evolution (jHollandl . Il975f ). Starting from a randomly 



generated population of P rule tables /„, the algorithm optimizes possible so- 
lutions by evaluating the fitness function 

(T„-T).7Vc ^ ^ ^"'^^'^ 

" J t=T^-T \ i=0 

Nc-1 \ 

+ E {i-'^^rw^^} . (32) 

i=[a-Nc] / 

The optimization algorithm then is defined as follows: 

1. Generate a random initial population J- = {/i, /p} of rule tables. 

2. Randomly assign system sizes A^™"' < N" < N™"-'-^ to all rule tables. 

3. For each rule table, generate a random initial state vector. 

4. Randomly mutate one entry of each rule table (generating a population 
T* of mutants). 

5. Iterate dynamics over time steps for J- and T* . 

6. Evaluate $(/„) and <&(/i^) for all rule tables < n < P, averaging over 
the past Tu — T update steps (with an additional penalty term if /„ does 
not converge to a fixed point). 

7. For each n, replace /„ with /*, if $(/„) < $(/*). 

8. Replace the least fit solution by a duplicate of the fittest one. 

9. Go back to step 2 and iterate. 

The outcome of this search algorithm is a set of rule tables, which then can be 
"translated" into (spatially coupled) Boolean networks or threshold networks 
with suitable thresholds. This yields a set of (minimal) dynamical networks 
which solve the pattern formation task by means of internal information pro- 
cessing. 



A. 2. Evolution of cellular automata 

The genetic algorithm sketched above is run with the following parameter 
choices: 15 < Nc < 150, i.e. during GA runs the system size is varied randomly 
between 15 and 150 cells, and the population size is set to P = 100. Fig. [2T] 
shows the fitness of the highest fitness mutant and the average fitness of the 
population as a function of the number of successive mutation steps during op- 
timization. A useful solution is found rather quickly (after about 200 updates), 
with further optimization observed during further 10000 generations. At time 
step 10000 mutations are turned off, thus now the average fitness of the estab- 
lished population under random initial conditions and random fiuctuations of 
the system size Nc is tested. The average fitness $ sa 0.98 indicates a surpris- 
ingly high robustness against fluctuations in the initial start pattern, indicating 
that the system is capable of de novo pattern formation. In the "fitness picture" , 
Fig. [22] confirms that the dynamically regulated domain size ratio a/(l — a) 
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Figure 21: Average fitness {^)(t) of the mutant population JT* and fitness of tlie liigliest 
fitness mutant $best (t) as a function of simulation time during tlie genetic algorithm 
run that lead to the high fitness solution used in this paper. At time step 10000 
mutations were turned off, in order to test the established population of optimized 
rule tables under different initial conditions (this corresponds to the sharp increase 
of at time step 10000). The evolved population of rule tables has an average 

fitness of about 0.98, independent from the initial conditions and system size Nc (in 
the tested range, i.e. 15 < Nc < 150). 



indeed is independent of system size (proportion regulation), there is only a 
weak decay of the fitness at small values of Nc ■ Interestingly, one also observes 
that for regeneration of Hy dra polyps from random cell aggregates a minimum 
number of cells is required (jTechnau et al 1 12000 ). The model suggests that this 
observation might be explained by the dynamics of an underlying pattern gen- 
erating mechanism, i.e. that there has to be a minimum diversity in the initial 
condition for successful de novo pattern formation. 



B. Statistical analysis of solutions 

An interesting question is how "difficult" it would be for an evolutionary pro- 
cess driven by random mutations and selection to find solutions for the pattern 
formation problem based on neighbor interactions between cells. As we showed 
above, the genetic algorithm finds the correct solution fast, however, this does 
not necessarily mean that biological evolution could access the same solution as 
fast. If there is only one, singular solution, evolution may never succeed find- 
ing it, as the genotype which already exists cannot be modified in an arbitrary 
way without possibly destroying function of the organism [developmental con- 
strains). To illustrate this point, we generated an ensemble of Ne ~ 80 different 
solutions with <i> > 0.96 and performed a statistical analysis of the rule table 
structure. 
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Figure 22: Average fitness of the highest fitness rule table as a function of the system 
size Nc- For system sizes Nc > 80 the fitness is almost constant at about 0.98. Notice 
that the decrease of the fitness for small Nc is an effect of the dynamics, not of the 
genetic algorithm implementation (all Nc in the range 15 < Nc < 150 were tested 
with equal probability), hence the dynamics of pattern formation may impose a lower 
boundary on the range of system (animal) sizes tolerated by natural selection. 

As one can see in Fig. 1231 some positions in the rule table are quite fixed, 

1. e., there is not much variety in the outputs, whereas other positions are more 
variable. We note that a number of rule table positions are a priori fixed due 
to the constraints imposed on dynamics. For example, to support a stable 
boundary ...2221000... as for the model described, five rules become fixed: 222 

2, 221 2, 210 1, 100 and 000 (these rules can be clearly 
distinguished as pronounced peaks in Fig. [23]). Consistency with boundary 
conditions fixes more rules. For the model as described in this manuscript, i.e. 
boundary conditions (t_i — (t^^ — 0, the rule 022 2 becomes fixed, too. 

However, the output frequency distribution alone does not allow to really 
judge the "evolvability" of the solutions: if there are strong correlations between 
most of the rule table entries, evolutionary transitions from one solution to 
another would be almost impossible. To check this point, we studied statistical 
two point correlations between the rule table entries. The probability for finding 
state a at position a and state a' at rule table position h is given by 



where 6 is the Kronecker symbol and n runs over the statistical ensemble of size 




11=1 



(33) 
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Figure 23: Frequency distribution p{(Ji) of outputs as a function of the rule table index 
as denoted in table 1. Ensemble statistics is taken over 80 different solutions with 
$ > 0.96. The upper panel shows the distribution for ai = 0, the middle panel the 
distribution for cr,; — 1 and the lower panel the distribution for cr,; — 2. 



Ne ■ The two point correlation between a and b then is defined as 

C"'' = ci f max p"'' (a, cr' ) - C2 ) (34) 

with ci = 9/8 and C2 = 1/9 to obtain a proper normalization with respect to 
the two limiting cases of equal probabilities (p'^^{a,a') = 1/9 V(cr, cr')) and 
p»f'(cr, cr') = 1 for cr = 5-, a = a and p"''(cr, cr') = for all other (cr, cr')). Fig. 
[24l shows the frequency distribution of C°''(cr, a'), averaged over all possible pairs 
(a, 6). About 65% of rule table positions are strongly correlated (C"'' = 1.0), the 
rest shows correlation values between 0.3 and 1.0. Hence, we find that the space 
of solutions is restricted, nevertheless there is variability in several rule table 
positions. To summarize this aspect, the pattern formation mechanism studied 
in this paper shows considerable robustness against rule mutations, however, 
a "core module" of rules is always fixed. Interestingly, a similar phenomenon 
is observed in developmental biology: Regulatory modules involved in develop- 
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Figure 24: Frequency distribution p{C°'') of two point correlations C°^{a,a ) of rule 
table entries, as defined in Eqn. (B2), averaged over all possible pairs of rule table 
entries. About 65% of rule table entries have correlation 1.0, for the rest the correlation 
is lower. 



mental processes often are ev olutionarily very conservative , i.e., they are shared 
by almost all animal phyla (jPavidsonl . 2001), while morphological variety is 
created by (few) taxon specific genes (jBoscbl . 120031 ) and rewiring of existing 
developmental modules. 



C. Numerical analysis of the dynamical transition at — 

The rate equation of the pattern formation system in the presence of noise 
with error rate r^, is given by 

2a*re = (1 - a*)re, (35) 

i.e. a* — 1/3. Comparison to the equilibrium position in the noiseless case 
indicates that the system undergoes a step-like discontinuity with respect to a* 
at Te = 0. A numerical analysis that considers small variations of Ve close to 
zero and averages over time windows of variable length Sw can be applied for 
supporting numerical evidence and finite size scaling. 

Figs . [25l and [26l show noise dependence and finite size scaling of the transition 
from the unperturbed solution to the solution under noise. Considering update 
time windows of different length s^, in case of a discontinuity (i.e., a step- 
like 'jump' of the order parameter) at re = we would expect a shift of the 
transition point r*'''*"''(sm) towards = which is proportional to as well 
as a divergence of the slope at the transition point when s^ is increased, i.e. 
da* /dre{rl^°'^'') oo when Sw ^ oo. 



37 




Figure 25: Average domain boundary position a* as a function of the error rate r^, 
sampled over update windows of different lengths (ensemble statistics, 400 different 
initial conditions for each data point). The abscissa is logarithmic. With increasing 
Sw, the transition from the solution q^^j = 0.281 under deterministic dynamics to 
Q* = 1/3 under noise is shifted towards =0. The two straight lines define a lower 
boundary Q-l^^ and a upper boundary Q*p, as explained in the text. 



The shift of r*''°"''(su,) is most easily measured by defining a lower and a 
upper boundary a^^^ and a*^, respectively (Fig. [^5]) : when a* crosses these 
boundaries, two transition points r"^ and r^"™ are obtained. We find that 

~ CupS~^ and r^°^ « ciowS~^ with Cup > ciow as expected (Fig. which 
imphes that the diflFerence Ar*''°"*(suj) :— r"P — r^"™ scales as 

Ar*™"^(s^) = {cup - cio^) s-\ (36) 

hence, because Aa*(r*''''"'') = const. = a*^ - a'l^.^^, indeed da*/c?re(r*''""*) 
diverges when the sampling window size goes to infinity. 



D. Derivation of stability conditions for growing systems 

D.l. Symmetric growth at the system boundaries 

Let us assume we start with a system of Nq cells, with an initial boundary 
position at cell Ni. In the deterministic case rg = 0, it is straight-forward to 
see that the time dependence of the average boundary position is given by 

a-m - (3T) 

Na+rgt 

hence we have 

= lim a*it) = i. (38) 
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Figure 26: Finite size scaling of the upper and lower transition points r"*" and r'f^, i.e. 
the points where a* crosses ai^^ and aZp, respectively (Fig. [25}, as a function of the 
sampling window length s^. Both r"*" and r'j^^ vanish oc s^^, as indicated by the line 
with slope —1 in this log- log-plot. 



In the case re > 0, the time dependence of the average boundary position is 
given by 

m + ^rgt~ 2a*{t)r,t + (1 - a*{t))r,t 
a (t) ~ — , (39) 



which simplifies to 



a*it) = + (40) 

^ ' No + rgt + 3ret ^ ' 



Assuming infinite growth, this leads to the asymptotic boundary position 



hma*(i) = %±^. (41) 

Fig. [TBI shows a°°{rg) for four different values of re; it becomes evident that an 
approximately 'correct' proportion regulation requires to be at the order of 
rg or larger, i.e. r^/rg > 1. While (the rate of regulatory signals) may not be 
increased significantly above the growth rate Vg, due to metabolic constraints, 
in later stages of development the steady decrease of rg will ensure that the 
condition r^/rg > 1 is fuUfiUcd and proportion regulation approaches the steady 
state of the adult organism. 

D.2. Growth by homogeneous cell proliferation 

We require that de-novo pattern formation has taken place in a system of 
stationary size A^o and has converged to its final pattern. 
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Assuming homogeneous cell proliferation with a probability pd per cell, it 
is easy to see that proliferation of "blue" {ai = 2) and "black" {at = 0) cells 
conserves pattern proportions, proliferation of the "red" boundary cell {ai = 1), 
however, leads to readjustment of the boundary at its original position before 
proliferation, and hence slightly reduces a (see inset of Fig fTT)) . Since, in this 
case, boundary readjustment needs two update time steps and occurs with prob- 
ability Pd, a{t) is given by 

N,it-l)il+pd)-pd/2 

= m • ^ ^ 

System size grows geometrically, i.e. N{t) ~ No{l +PdY- Inserting this de- 
pendence into Eq. dH]) and using Ni{t - 1) = a{t ~ l)N{t - 1), it follows 
that 

Recursively inserting for a{t — r) for r e {1, t — 1}, we conclude that a{t) is 
given by 

«W = "0 - H TS^+P'^y^ ^ "° " 2^ " i^+Pdrl (44) 

This implies that the asymptotic boundary position in the limit t ^ cx) is 
independent from pd'. 

Hence deviations at the order of I/Nq from the boundary position ao = 0.281 
of non-growing systems are found. 

For the case when noise is present, it is not possible to find a general solution 
since proliferation can affect both the velocity and the type of particles travelling 
through the domains in intricate ways (in the case of symmetric growth at the 
system boundaries, as discussed in the previous subsection, this problem is 
avoided). If is very small, however, we can assume that proliferation events 
and particle propagation are essentially decoupled and that the system has 
enough time to relax to a stationary state between proliferation events 0. In 
this limit, we can generalize Eqn. psp in a straight-forward way: 

2are + Pd = {I - a)re, (46) 

leading to 



•^Furthermore, one can show that single proUfcration events that occur in or near quasi- 
particles have the same effect as certain subclasses of one-site errors as discussed in section 
|4] From this we conclude that - for moderate - the statistics of boundary readjustments 
remains unchanged, just as it was found in section |4] for error rates re < 1/2. 
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Te = pN is a monotonously growing function in time for fixed error probability 
p, hence it follows that, for large t, the boundary position converges to the same 
value a = 1/3 as for constant size systems. 
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